Informe autocontenido Análisis Exploratorio de Datos Monitoreo de la pesquería Bacalao de profundidad APA

Analisis exploratorio 1986-2021

Author

Mauricio Mardones I

Published

July 20, 2023

ANTECEDENTES

El presente reporte tiene un codigo autocontenido con los Analísis Exploratorios de Datos (AED) de la pesquería de bacalao de profundidad extraído por la flota pesquera artesanal (espinel). Los analisis aca descritos estan compuestos de dos bases. La primera es de los registros oficiales de captura recogidos por el SERNAPESCA. En segundo lugar, se presentan el AED de los datos del monitoreo d la pesqería llevado a cabo por el IFOP a través del Programa de Seguimiento de Pesquerías Demersales del Departamento de Evaluación de Pesquerías.

El principal objetivo es identificar vacios y fortalezas de los datos en terminos espacio-temporales, y con ello tomar decisiones para el paso posterior de modelación de la dinámica poblacional del recurso para realizar asesoría en terminos de estatus y recomendación de una Captura Biologicamente Aceptable (CBA).

ANÁLISIS EXPLORATORIO DE DATOS (AED)

Cargo librerías necesarias para el análisis exploratorio de los datos de las distintas bases de bitácora, tallas y biológico.

library(here)
#analisis
library(ggsignif)
library(ggrepel)
library(ggpubr)
library(inlmisc)
library(nortest) #para testear distribucion
library(skimr) #provides a frictionless approach to summary statistics w
library(easystats) # multiples unciones analiticas
library(lme4)
library(skimr)
library(readxl)
library(fitdistrplus)
# vizualizacion
library(ggridges)
library(sf)
library(GGally)
library(tidyverse, quietly = TRUE)
library(knitr, quietly = TRUE)
library(kableExtra)
library(raster)
library(egg)
library(car) #Variance inflation Factor
library(ggthemes)
library(sjPlot)
library(GGally)
# A function for dotplots
multi_dotplot <- function(filename, Xvar, Yvar){
  filename %>%
    ggplot(aes(x = {{Xvar}})) +
    geom_point(aes(y = {{Yvar}}),
               alpha=0.4) +
    theme_bw() +
    coord_flip() +
    labs(x = "Order of Data")}

Identifico los directorio de trabajo y leo las tres bases de datos; a saber:

  • Datos FIP (Rendimiento) 96-32
  • Bitácora
  • Biologico
  • Tallas
  • Desembarques (1985-2022)
bit <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BITACORAS ESPINEL bacalao_corr.xlsx", sheet = "BITACORAS_ESPINEL_bacalao",  guess_max = 100000)
long <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/LONGITUD_ESPINEL_bacalao.xlsx")
bio <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BIOLOGICO ESPINEL bacalao.xlsx", 
    sheet = "BIOLOGICO_ESPINEL_bacalao")
landing <- read_excel("Artesanal_Historica/DESEMBARQUE HISTÓRICO.xlsx", 
    skip = 2)
bit8696<- read_csv2("bit_art_86_20.csv")

Identificamos los registros asociados recurso objetivo

  • Codigo del Recurso: 37
  • Nombre común: Bacalao de profundidad
  • Nombre científico: Dissostichus eleginoides

Idetifico los registros por recurso y por base para luego filtrar.

dim(bit)
   [1] 14208    77
table(bit$COD_ESPECIE)
   
       2     4     5     6    14    15    16    17    18    20    22    23    24 
       7     1    14    25   147     6    89    16   698     1    25     3     7 
      25    27    33    35    37    42    43    49    53    56    59    74    75 
       3    11     2     2 11074     3     9    21     1     5    35     1    28 
      81    88    99   103   104   105   106   108   119   127   129   144   194 
      19     1   159    21    24     2     7    44     1     1     1     1    91 
     199   212   226   239   243   323   324   325   334   844   849   920   984 
     178    18     1    22     1     2    50     6   830     3   138     1     1 
     989   999  1309 
      36   249    58
dim(long)
   [1] 15831    24
table(long$COD_ESPECIE)
   
      37 
   15831
dim(bio)
   [1] 127409     36
table(bio$COD_ESPECIE)
   
       37 
   127409
# las bases de biologicos y longitud solo teienen  registros de bacalao.

Comenzar a trabajar bases por separado

DESEMBARQUES

Tabla con los desembarques oficiales (Sernapesca, 2022)


kbl(landing, booktabs = T,format = "html",
    caption = "Desembarque Bacalao Artesanal por Región") %>%
    kable_styling(latex_options = c("striped",
                                  "condensed","scale_down"),
                full_width = FALSE)
Desembarque Bacalao Artesanal por Región
...1 AyP TPCA ANTOF ATCMA COQ VALPO LGBO MAULE ÑUBLE BBIO ARAUC RIOS LAGOS AYSEN MAG ...17
1985 NA 8 102 146 162 1674 NA 349 NA 1527 NA NA 63 28 NA 4059
1986 NA 842 1531 213 329 1242 NA 88 NA 1546 NA NA 311 6 NA 6108
1987 NA 182 335 171 51 671 NA NA NA 1264 NA NA 689 21 NA 3384
1988 NA 131 174 251 42 670 NA 46 NA 1597 NA NA 877 8 NA 3796
1989 NA 217 175 177 109 756 NA 114 NA 1859 10 NA 1462 8 NA 4887
1990 NA 335 161 126 100 444 NA 5 NA 2620 NA NA 1789 36 NA 5616
1991 NA 280 186 142 223 519 NA NA NA 1901 NA NA 680 NA NA 3931
1992 NA 30 45 77 53 477 NA 131 NA 2108 1 NA 741 1 NA 3664
1993 NA 37 4 66 81 1676 NA 295 NA 1111 NA NA 807 NA 45 4122
1994 NA 24 139 213 68 671 NA 526 NA 1920 NA NA 1374 41 411 5387
1995 NA 252 259 286 65 206 NA 190 NA 1347 NA NA 1680 31 266 4582
1996 NA 476 220 121 177 272 NA 226 NA 1390 2 NA 1989 22 83 4978
1997 NA 36 82 47 57 98 NA 116 NA 1077 1 NA 1803 22 83 3422
1998 NA 9 67 169 61 155 NA 399 NA 1284 NA NA 2038 NA 11 4193
1999 NA 42 72 168 87 484 NA 406 NA 1369 NA NA 2910 NA 270 5808
2000 NA 285 98 153 105 200 NA 292 NA 800 NA NA 3515 NA 345 5793
2001 NA 198 85 127 113 191 NA 136 NA 1049 NA NA 2045 NA NA 3944
2002 NA 154 44 97 57 188 NA 186 NA 732 NA NA 3097 NA 10 4565
2003 NA 58 25 51 27 338 NA 123 NA 752 NA NA 3189 NA 179 4742
2004 NA 76 32 63 23 308 NA 132 NA 412 NA NA 2303 NA 70 3419
2005 NA 111 14 47 26 84 NA 175 NA 475 NA NA 2021 NA 325 3278
2006 NA 87 3 50 50 92 NA 159 NA 243 NA NA 1367 NA 40 2091
2007 19 43 7 22 25 43 NA 84 NA 189 NA 501 1157 NA NA 2090
2008 4 48 21 9 3 32 NA 70 NA 219 NA 349 803 NA NA 1558
2009 13 85 24 20 16 79 NA 112 NA 239 NA 522 571 NA NA 1681
2010 13 43 13 58 11 37 NA 63 NA 109 NA 465 655 NA NA 1467
2011 19 82 20 46 6 49 NA 90 NA 276 NA 613 987 NA NA 2188
2012 25 68 11 82 14 26 NA 105 NA 317 NA 290 1126 NA NA 2064
2013 26 58 12 68 23 52 NA 214 NA 186 NA 277 622 NA 20 1558
2014 13 48 - 59 34 64 NA 171 NA 129 NA 252 453 NA 57 1280
2015 6 52 13 53 38 93 NA 119 NA 261 NA 298 638 NA 38 1609
2016 17 51 20 59 16 70 NA 71 NA 335 NA 327 755 NA NA 1721
2017 20 92 33 89 33 117 NA 116 NA 316 NA 312 790 NA 28 1946
2018 29 78 5 74 49 78 NA 113 NA 251 NA 218 620 NA 25 1540
2019 26 55 7 63 35 42 NA 102 NA 248 NA 203 696 NA 94 1571
2020 5 18 7 NA 1 27 NA 78 NA 145 NA 229 261 NA 32 803
2021 2 83 10 138 17 52 NA 206 NA 304 NA 500 517 NA 21 1850
2022 26 100 - 148 58 75 - 143 - 318 - 332 754 - 29 1983

Ploteo los desembarques por region y por año

landing <- as.data.frame(lapply(landing, as.double))
str(landing)
   'data.frame':    38 obs. of  17 variables:
    $ ...1 : num  1985 1986 1987 1988 1989 ...
    $ AyP  : num  NA NA NA NA NA NA NA NA NA NA ...
    $ TPCA : num  8 842 182 131 217 335 280 30 37 24 ...
    $ ANTOF: num  102 1531 335 174 175 ...
    $ ATCMA: num  146 213 171 251 177 126 142 77 66 213 ...
    $ COQ  : num  162 329 51 42 109 100 223 53 81 68 ...
    $ VALPO: num  1674 1242 671 670 756 ...
    $ LGBO : num  NA NA NA NA NA NA NA NA NA NA ...
    $ MAULE: num  349 88 NA 46 114 5 NA 131 295 526 ...
    $ ÑUBLE: num  NA NA NA NA NA NA NA NA NA NA ...
    $ BBIO : num  1527 1546 1264 1597 1859 ...
    $ ARAUC: num  NA NA NA NA 10 NA NA 1 NA NA ...
    $ RIOS : num  NA NA NA NA NA NA NA NA NA NA ...
    $ LAGOS: num  63 311 689 877 1462 ...
    $ AYSEN: num  28 6 21 8 8 36 NA 1 NA 41 ...
    $ MAG  : num  NA NA NA NA NA NA NA NA 45 411 ...
    $ ...17: num  4059 6108 3384 3796 4887 ...
summary(landing)
         ...1           AyP             TPCA            ANTOF        
    Min.   :1985   Min.   : 2.00   Min.   :  8.00   Min.   :   3.00  
    1st Qu.:1994   1st Qu.:11.25   1st Qu.: 44.25   1st Qu.:  12.75  
    Median :2004   Median :18.00   Median : 77.00   Median :  32.50  
    Mean   :2004   Mean   :16.44   Mean   :128.26   Mean   : 112.67  
    3rd Qu.:2013   3rd Qu.:25.25   3rd Qu.:148.25   3rd Qu.: 111.25  
    Max.   :2022   Max.   :29.00   Max.   :842.00   Max.   :1531.00  
                   NA's   :22                       NA's   :2        
        ATCMA            COQ             VALPO             LGBO    
    Min.   :  9.0   Min.   :  1.00   Min.   :  26.0   Min.   : NA  
    1st Qu.: 58.0   1st Qu.: 23.50   1st Qu.:  65.5   1st Qu.: NA  
    Median : 82.0   Median : 49.50   Median : 136.0   Median : NA  
    Mean   :106.7   Mean   : 64.34   Mean   : 325.1   Mean   :NaN  
    3rd Qu.:148.0   3rd Qu.: 77.75   3rd Qu.: 468.8   3rd Qu.: NA  
    Max.   :286.0   Max.   :329.00   Max.   :1676.0   Max.   : NA  
    NA's   :1                                         NA's   :38   
        MAULE           ÑUBLE          BBIO            ARAUC           RIOS      
    Min.   :  5.0   Min.   : NA   Min.   : 109.0   Min.   : 1.0   Min.   :203.0  
    1st Qu.: 99.0   1st Qu.: NA   1st Qu.: 253.5   1st Qu.: 1.0   1st Qu.:270.8  
    Median :127.0   Median : NA   Median : 603.5   Median : 1.5   Median :319.5  
    Mean   :165.3   Mean   :NaN   Mean   : 848.0   Mean   : 3.5   Mean   :355.5  
    3rd Qu.:194.0   3rd Qu.: NA   3rd Qu.:1363.5   3rd Qu.: 4.0   3rd Qu.:473.8  
    Max.   :526.0   Max.   : NA   Max.   :2620.0   Max.   :10.0   Max.   :613.0  
    NA's   :2       NA's   :38                     NA's   :34     NA's   :22     
        LAGOS            AYSEN            MAG             ...17     
    Min.   :  63.0   Min.   : 1.00   Min.   : 10.00   Min.   : 803  
    1st Qu.: 661.2   1st Qu.: 8.00   1st Qu.: 28.25   1st Qu.:1753  
    Median : 842.0   Median :22.00   Median : 51.00   Median :3402  
    Mean   :1267.2   Mean   :20.36   Mean   :112.82   Mean   :3228  
    3rd Qu.:1799.5   3rd Qu.:29.50   3rd Qu.:157.75   3rd Qu.:4472  
    Max.   :3515.0   Max.   :41.00   Max.   :411.00   Max.   :6108  
                     NA's   :27      NA's   :16
landing2 <- landing %>% 
  pivot_longer(cols=c("AyP"  , "TPCA"  ,"ANTOF", "ATCMA", "COQ",
               "VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
               "RIOS" , "LAGOS", "AYSEN", "MAG"), 
               names_to = "REGION", 
               values_to = "CAPTURA") %>% 
  rename(AÑO="...1") %>% 
  dplyr::select(-2)

Genero un gráfico de barras


landing2$REGION <- factor(landing2$REGION , 
                          levels = c("AyP"  , "TPCA"  ,"ANTOF", "ATCMA", "COQ",
               "VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
               "RIOS" , "LAGOS", "AYSEN", "MAG"))
desem <- ggplot(landing2 %>% 
                  filter(REGION != "MAG") %>% 
                   drop_na(REGION),aes(AÑO, CAPTURA, fill=REGION)) +
  geom_bar(stat="identity")+
  scale_fill_viridis_d(option="G")+
  theme_few()+
  scale_x_continuous(breaks = seq(from = 1985, to = 2022, by = 5))+
  scale_y_continuous(breaks = seq(from = 0, to = 4000, by = 1000))+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank(),
        legend.position = "none")+
  facet_wrap(~REGION, ncol=7)+
  labs(y="Desembarques Oficiales (t)",
       x="")
desem

Los principales desembarques estan asociados a las regiones de Antifagasta, Valparaiso y BioBIo pero solo durante los primeros años (1985- mediados del 2000). Luego de esto, todas las regiones vieron disminuidos los registros de extracción del recurso.

BITÁCORA

Esta base tiene como principal objetivo obtener un indicador de esfuerzo para calcular un indice de abundancia relativo como la CPUE.

Primero identifico la estructura de la base y filtro el rrecurso bacalao

# filtro bitacoras dejando solo bacalao
bitb <- bit %>% 
  filter(COD_ESPECIE==37)%>% 
  mutate(dfp = as.double(dfp), 
         PESO_corr = as.double(PESO_corr),
         prof_med=as.double(prof_med))
dim(bitb)
   [1] 11074    77

se eliminan 1208-11074 = 3134 registros

glimpse(bitb)

ahora las estadísticas descriptivas de las variables de intéres. En este caso dfp, PESO_corr, prof_med

Miramos los oultiers de los datos

# OUTLIERS

#Order data
bitb <- bitb %>%
  mutate(order = seq(1:nrow(bitb)))

#Select continuous variables to plot
p1 <- multi_dotplot(bitb, order, dfp)
p2 <- multi_dotplot(bitb, order, PESO_corr)
p3 <- multi_dotplot(bitb, order, prof_med)
p4 <- multi_dotplot(bitb, order, N_TRIPULANTES)
p5 <- multi_dotplot(bitb, order, HORAS_REPOSO)
p6 <- multi_dotplot(bitb, order, NUMERO_DE_ANZUELOS)

#Plot as a grid
ggarrange(p1, p2, p3, p4,p5, p6,  nrow = 3)

Remuevo outliers

Los datos de profundidades andan bien, Filtro los datos de la variable dfp entre 0 y 30 días. y repito la inspección

bitb2 <- bitb  %>% 
  filter(dfp>0,
         dfp<80,
         PESO_corr<50000,
         N_TRIPULANTES<30,
         NUMERO_DE_ANZUELOS<40000)

Miro nuenvamente los datos filtrados

# OUTLIERS

#Order data
bitb2 <- bitb2 %>%
  mutate(order = seq(1:nrow(bitb2)))

#Select continuous variables to plot
p1a <- multi_dotplot(bitb2, order, dfp)
p2a <- multi_dotplot(bitb2, order, PESO_corr)
p3a <- multi_dotplot(bitb2, order, prof_med)
p4a <- multi_dotplot(bitb2, order, N_TRIPULANTES)
p5a <- multi_dotplot(bitb2, order, HORAS_REPOSO)
p6a <- multi_dotplot(bitb2, order, NUMERO_DE_ANZUELOS)

#Plot as a grid
ggarrange(p1a, p2a, p3a, p4a, p5a, p6a,  nrow = 3)

Identifico las dimensionesd e la nueva base y la comparo con la anterior

dim(bitb2)
   [1] 4745   78
dim(bitb)
   [1] 11074    78
table(bitb$año)
   
   1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 
     15  557 1046 1072  585  931  576  196  131   83   43   31   57   47   73   20 
   2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 
     85  497  506  583  760 1175  474  497  564  470

Los registros se minimiza en ciertos años lo cual genera problemas en la fiabilidad de la estimación.

ahora calculamos el rendimiento nominal

CPUE Nominal

Ahora estimo la CPUE, la cual es el rendimoento entre el esfuerzo y la captura. La medida de esfuerzo será dfp. También calculo un rendimiento con la variable de esfuerzo HORAS_REPOSO pero solo como antecedente dado que no tengo muchos datos.

bitb2 <- bitb2 %>% 
  mutate(CPUE = PESO_corr/dfp) %>% 
  mutate(CPUE2 =PESO_corr/HORAS_REPOSO)
bitb2$prof_med_cat <- cut(bitb2$prof_med, 
                          breaks = c(500, 1000, 1500,
                                     2000, 3000))
# Frequency polygon plot for catch
his <- ggplot(bitb2, aes(CPUE)) +
  geom_freqpoly(bins = 5) +
  labs(x = "dfp", y = "Frequency") +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", 
                                    fill=NA, size = 1))+
  facet_wrap(~año)+
  theme_few()
# Patterns in the variance? (any lack of homogeneity)
point <- ggplot(bitb2 %>% 
                  filter(PESO_corr>1000), aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 5, alpha = 0.6) +
  stat_smooth(method = "lm")+
  theme(panel.background = element_blank()) +
  theme(panel.border = element_rect(fill = NA, size = 1)) +
  theme(strip.background = element_rect(fill = "white", 
                   color = "white", size = 1)) +
  theme(text = element_text(size=13)) +
  xlab("dfp") + ylab("Peso Corregido")


#Plot as a grid
ggarrange(his, point,  nrow = 1)

Tendencia del esfuerzo dfp a través de los años y por región

effor <- ggplot(bitb2 %>%
                  group_by(año, REGION_PUERTO_RECALADA) %>% 
                  summarise(DFPM=mean(CPUE)),
                aes(año, DFPM))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 3))+
  geom_smooth(method = 'lm', 
              colour = 'blue', 
              size = 1.5)+
  facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="Effort (dfp)",
       x="")+
  ylim(5,500)
effor

efformean <- bitb %>% 
  group_by(año) %>% 
  summarise(dfpmean= mean(dfp))
write.table(efformean, "dfptable.txt")

Tendencia del esfuerzo NUMERO_DE_ANZUELOS a través de los años y por región

efforanz <- ggplot(bitb2 %>%
                  group_by(año, REGION_PUERTO_RECALADA) %>% 
                  summarise(NANZ=mean(NUMERO_DE_ANZUELOS)),
                aes(año, NANZ))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
  geom_smooth(method = 'lm', 
              colour = 'green', 
              size = 1.5)+
  facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="Effort (dfp)",
       x="")
efforanz

Calculo los ceros de dfp

#CALCULATE NUMBER OF ZEROS

# What is the percentage of zeros i the response variable

round(sum(bitb2$dfp == 0) * 100 / nrow(bitb2),0)
   [1] 0
#0

Interacciones

# Interactions

# Year x season
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("Points") + ylab("Catch") +
  facet_grid(año~trim)

# No

ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~N_TRIPULANTES)+
  theme_few()

# Perhaps

ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  theme_bw() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(trim~N_TRIPULANTES)+
  theme_few()

# CPUE slope varies between habitats - interaction

Profundidad por region

profu <- ggplot(bitb2, aes(año, desc(prof_med), 
                            group=año))+
  geom_boxplot(fill=NA, alpha=.5)+
  geom_jitter(size=0.4, alpha=0.2,
              width = .25)+
  facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
  theme_few()+
  scale_x_continuous(breaks = seq(from = 2011, to = 2022, by = 3))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  labs(x="Región",
       y="Profundidad (mts)")
profu

Comparo las medias de la profundidad

# Prueba de Kruskal-Wallis
pairwise.t.test(x = bitb2$prof_med, 
                g = bitb2$REGION_PUERTO_RECALADA, 
                p.adjust.method = "holm",
                pool.sd = TRUE, 
                paired = FALSE, 
                alternative = "two.sided")
   
    Pairwise comparisons using t tests with pooled SD 
   
   data:  bitb2$prof_med and bitb2$REGION_PUERTO_RECALADA 
   
      1       3       4       5       7       8       10      12      14     
   3  1.00000 -       -       -       -       -       -       -       -      
   4  1.00000 1.00000 -       -       -       -       -       -       -      
   5  3.2e-07 0.09875 1.00000 -       -       -       -       -       -      
   7  1.00000 1.00000 1.00000 8.0e-11 -       -       -       -       -      
   8  5.5e-08 0.02686 1.00000 1.00000 5.2e-11 -       -       -       -      
   10 0.00032 0.43592 1.00000 1.00000 0.00012 1.00000 -       -       -      
   12 0.68271 0.27691 0.21334 0.01917 0.44715 0.01135 0.02686 -       -      
   14 1.00000 1.00000 1.00000 0.00156 1.00000 0.00027 0.04753 0.35077 -      
   15 2.1e-12 1.2e-13 1.7e-07 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000 9.3e-15
   
   P value adjustment method: holm
ggplot(bitb2, 
       aes(x = nombre_puerto_recalada , 
           y = dfp, 
           fill = nombre_puerto_recalada )) +
  geom_violindot(fill_dots = "black") +
  theme_modern() +
  scale_fill_material_d()


ggplot(bitb2, aes(x = as.factor(REGION_PUERTO_RECALADA ), 
                  y = bitb2$NUMERO_DE_ANZUELOS, 
                  fill = as.factor(REGION_PUERTO_RECALADA))) +
  geom_violin(scale="width") +
  theme_modern() +
  scale_fill_material_d(palette="ice",
                        name="REGION")

ggplot(bitb2 %>% 
         drop_na(prof_med_cat), 
       aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', 
              colour = 'red', 
              se = FALSE) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(prof_med_cat~año)

ploteo los datos de CPUE totales.

cpuenom <- ggplot(bitb2 %>% 
  group_by(zona, año) %>%
  summarise(CPUEM=mean(CPUE)) %>% 
  filter(CPUEM<600,
         año>1996,
         zona!=4),
  aes(año, CPUEM, color=zona, group=zona))+
    geom_point()+
  geom_smooth()+
  scale_color_viridis_b(option="B")+
  scale_x_continuous(breaks = seq(from = 2010, to = 2022, by = 1))+
  theme_few()+
  labs(x="",
       y="CPUE")
cpuenom

guardo los datos bitacora CPUE nominal 1996-2022 del objeto cpuenom

write_csv(cpuenom, "CPUE_MON.csv")

CPUE Estandarizada

Identificar los proncipales factores para modelar la variable.

Luego de una reunión con Seguimiento (Patricio Galvez) se indicaron algunos aspectos que deben ser considerados en la estandariación.

  • El poder de pesca (tipo de embarcación) es mayor en la zona centro sur de Chile. Identificar equilibrio de la base bitb2 respecto a este dato.

  • Usar zonas como factor en desmedro de región.

Identifico una estadistica descriptiva general de la base de estandarización

Correlaciones

Análisis de utilidad para identificar a traves de un metodo de corrrelación de pearson, la correlación en tre las variables que serán utilizadas en la estandariacion de la CPUE.

Primero identifico las variables a correlacionar

names(bitb2)
    [1] "ID"                           "embarque"                    
    [3] "COD_BARCO"                    "FECHA_HORA_RECALADA"         
    [5] "FECHA_HORA_ZARPE"             "COD_PESQUERIA"               
    [7] "COD_PESQUERIA_ant"            "N_TRIPULANTES"               
    [9] "N_CALADAS"                    "PUERTO_ZARPE"                
   [11] "PUERTO_RECALADA"              "nombre_puerto_recalada"      
   [13] "NRO_FORMULARIO"               "REGION_PUERTO_RECALADA"      
   [15] "ESPECIE_OBJETIVO"             "ECOSONDA"                    
   [17] "GPS"                          "POTENCIA_MOTOR_EQ"           
   [19] "VIRADOR"                      "GASTOS_VIAJE"                
   [21] "NUMERO_LANCE_EX"              "FECHA_LANCE"                 
   [23] "año"                          "mes"                         
   [25] "trim"                         "dfp"                         
   [27] "ID_CUADRICULA"                "ID_PROCEDENCIA"              
   [29] "PESO_TOTAL_CAPTURA"           "LATITUD"                     
   [31] "LONGITUD"                     "lat_ctgr"                    
   [33] "lon_ctgr"                     "zona"                        
   [35] "ESPECIE_OBJETIVO_LANCE"       "CLASE_LANCE"                 
   [37] "COD_ESPECIE"                  "PESO"                        
   [39] "PESO_corr"                    "CAPTURA_RETENIDA"            
   [41] "VOLUMEN_DESCARTE"             "LUGAR_DESCARTE"              
   [43] "TIPO_DESCARTE"                "PRECIO_UNITARIO"             
   [45] "DESTINO_CAPTURA"              "NRO_ANZUELOS"                
   [47] "REGION_PROCEDENCIA"           "FECHA_HORA_FIN_CALADO"       
   [49] "FECHA_HORA_INI_CALADO"        "LATITUD_INICIO_CALADO"       
   [51] "LONGITUD_INICIO_CALADO"       "LATITUD_FIN_CALADO"          
   [53] "LONGITUD_FIN_CALADO"          "FECHA_HORA_INI_VIRADO"       
   [55] "FECHA_HORA_FIN_VIRADO"        "LATITUD_INICIAL_VIRADO"      
   [57] "LONGITUD_INICIAL_VIRADO"      "LATITUD_FIN_VIRADO"          
   [59] "LONGITUD_FIN_VIRADO"          "PROFUNDIDAD_MINIMA_LM"       
   [61] "PROFUNDIDAD_MAXIMA_LM"        "PROFUNDIDAD_FONDO_VIRADO_INI"
   [63] "PROFUNDIDAD_FONDO_VIRADO_FIN" "prof_med"                    
   [65] "SEPARACION_ANZUELO"           "MARCA_ANZUELO"               
   [67] "NUMERO_DE_ANZUELOS"           "PERDIDA_ANZUELOS"            
   [69] "TIPO_CARNADA"                 "TAMANIO_ANZUELO"             
   [71] "INTENSIDAD_VIENTO"            "LONGITUD_LINEA_MADRE"        
   [73] "HORAS_REPOSO"                 "TIPO_ESPINEL"                
   [75] "N_PANOS"                      "ANZ_POR_PANO"                
   [77] "filtro"                       "order"                       
   [79] "CPUE"                         "CPUE2"                       
   [81] "prof_med_cat"
bitb3 <- bitb2 %>%
  dplyr::select(8, 18, 23,24, 25, 26, 32, 33, 34, 39, 69, 73, 75, 78, 79, 80, 81)
names(bitb3)
    [1] "N_TRIPULANTES"     "POTENCIA_MOTOR_EQ" "año"              
    [4] "mes"               "trim"              "dfp"              
    [7] "lat_ctgr"          "lon_ctgr"          "zona"             
   [10] "PESO_corr"         "TIPO_CARNADA"      "HORAS_REPOSO"     
   [13] "N_PANOS"           "order"             "CPUE"             
   [16] "CPUE2"             "prof_med_cat"
results <- correlation(bitb3)

results %>%
  summary(redundant=TRUE) %>% 
  plot(results, show_data = "points")+
  theme_bw()

skimr::skim(bitb3)
Data summary
Name bitb3
Number of rows 4745
Number of columns 17
_______________________
Column type frequency:
character 1
factor 1
numeric 15
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
TIPO_CARNADA 122 0.97 1 5 0 8 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
prof_med_cat 939 0.8 FALSE 4 (1e: 1867, (50: 961, (1.: 853, (2e: 125

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
N_TRIPULANTES 0 1.00 7.23 1.62 2.00 6.00 7.00 9.00 10.00 ▁▂▆▇▆
POTENCIA_MOTOR_EQ 167 0.96 405.74 399.03 0.00 320.00 390.00 420.00 3903.00 ▇▁▁▁▁
año 0 1.00 2018.05 2.47 2011.00 2016.00 2018.00 2020.00 2022.00 ▁▅▆▆▇
mes 0 1.00 7.29 3.79 1.00 3.00 9.00 11.00 12.00 ▅▂▁▂▇
trim 0 1.00 2.71 1.32 1.00 1.00 3.00 4.00 4.00 ▆▂▁▂▇
dfp 0 1.00 16.73 6.62 0.67 13.02 16.79 19.79 50.67 ▂▇▂▁▁
lat_ctgr 51 0.99 -35.49 6.89 -55.17 -40.99 -35.25 -33.17 -18.33 ▁▅▇▂▂
lon_ctgr 51 0.99 -73.24 1.62 -81.89 -74.69 -73.08 -71.97 -51.94 ▁▇▁▁▁
zona 0 1.00 2.11 0.65 1.00 2.00 2.00 2.00 4.00 ▂▇▁▃▁
PESO_corr 0 1.00 755.69 1304.56 0.00 41.00 110.00 1079.66 14092.00 ▇▁▁▁▁
HORAS_REPOSO 2988 0.37 1629.72 832.08 12.00 1200.00 1200.00 2400.00 7200.00 ▇▅▁▁▁
N_PANOS 1833 0.61 226.33 280.51 1.00 1.00 2.00 420.00 900.00 ▇▁▃▁▁
order 0 1.00 2373.00 1369.91 1.00 1187.00 2373.00 3559.00 4745.00 ▇▇▇▇▇
CPUE 0 1.00 54.95 92.65 0.00 2.44 6.49 84.01 1067.68 ▇▁▁▁▁
CPUE2 2988 0.37 1.55 5.23 0.00 0.51 0.99 1.79 158.44 ▇▁▁▁▁

Datos 0

Are there missing values? y sirve para identificar los factores que se utilzaran en la estandarización.

dim(bitb3)
   [1] 4745   17
colSums(is.na(bitb3))
       N_TRIPULANTES POTENCIA_MOTOR_EQ               año               mes 
                   0               167                 0                 0 
                trim               dfp          lat_ctgr          lon_ctgr 
                   0                 0                51                51 
                zona         PESO_corr      TIPO_CARNADA      HORAS_REPOSO 
                   0                 0               122              2988 
             N_PANOS             order              CPUE             CPUE2 
                1833                 0                 0              2988 
        prof_med_cat 
                 939

Es probable que HORAS_REPOSOde la red no podamos usarlo como factor ni como variable de esfuerzo por lo exiguo del dato.

Variabilidad

asd1<-ggplot(data = bitb3 %>% 
               filter(zona!=4), 
             aes(y=dfp, x=zona,
                 group=zona)) +
  geom_boxplot()+
  labs(title = "Variabilidad en DFP",
       y = "DFP", x = "ZONA") + 
  facet_wrap(~ año, scales = "free")+
   theme_few()
asd1

asd2<-ggplot(data =bitb3, 
             aes(x=LAT_CAT,
                 y=CPUE)) +
  geom_boxplot()+
  labs(title = "Variabilidad en distribución latitudinal",
       y = "CPUE", x = "Zona") + 
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  facet_wrap(~ año)
asd2

Normality and Homogenity factors

# Frequency polygon plot for catch
ggplot(bitb3 %>% 
         filter(zona!=4),
         aes(CPUE)) +
  geom_freqpoly(bins = 8) +
  labs(x = "Cpue", y = "Frequency") +
  facet_wrap(.~zona)+
  theme_few()+
  theme(panel.border = 
          element_rect(colour = "black", 
                       fill=NA, size = 1))

Shapiro-Wilk test for deviation from normality


shapiro.test(bitb3$CPUE)
   
    Shapiro-Wilk normality test
   
   data:  bitb3$CPUE
   W = 0.63039, p-value < 2.2e-16

Departure from normality…but not critical

Colinealidad entre factores

# COLLINEARITY
Coll <- c("año",  "dfp" ,"lat_ctgr", "zona" ,
          "N_PANOS" ,    "CPUE", "prof_med_cat")

# Obtain summary using the ggpairs command from the GGally library
ggpairs(bitb3[,Coll], ggplot2::aes(alpha = 0.9, colour=zona),
        lower = list(combo = "count"))+
  scale_fill_manual(values=c("red", "blue", "green",  "black")) +
   scale_colour_manual(values=c("red", "blue", "green", "black")) +
  theme_few()
# Nothing serious
#Calculate Variance Inflation Factor (VIF)
vifmodel <- round(vif(lm(PESO_corr ~ zona + dfp + año + mes,
                     data = bitb3)),2)
vifmodel
   zona  dfp  año  mes 
   1.34 1.29 1.06 1.01
barplot(vifmodel, main = "VIF Values", horiz = TRUE, col = "steelblue")

La multicolinealidad en el análisis de regresión se produce cuando dos o más variables predictoras están altamente correlacionadas entre sí, de modo que no proporcionan información única o independiente en el modelo de regresión.

Si el grado de correlación es lo suficientemente alto entre las variables, puede causar problemas al ajustar e interpretar el modelo de regresión.

La forma más común de detectar la multicolinealidad es mediante el factor de inflación de varianza (VIF), que mide la correlación y la fuerza de la correlación entre las variables predictoras en un modelo de regresión.

El valor de VIF comienza en 1 y no tiene límite superior. Una regla general para interpretar VIF es la siguiente:

Un valor de 1 indica que no hay correlación entre una variable predictora dada y cualquier otra variable predictora en el modelo. Un valor entre 1 y 5 indica una correlación moderada entre una variable predictora dada y otras variables predictoras en el modelo, pero esto a menudo no es lo suficientemente grave como para requerir atención. Un valor superior a 5 indica una correlación potencialmente grave entre una variable predictora dada y otras variables predictoras en el modelo. En este caso, es probable que las estimaciones de coeficientes y los valores p en el resultado de la regresión no sean fiables.

Interaccion

# Patterns in the variance? (Evidence for lack of homogeneity)
ggplot(bitb3, 
       aes(x = dfp, 
                  y = (CPUE))) +
  geom_jitter(shape = 16, 
              size = 2.5, 
              alpha = 0.3, 
              height = 0.25, 
              width = 0.5) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  facet_grid(bitb3$N_TRIPULANTES~bitb3$zona)+
  theme(panel.background = element_blank()) +
  theme(panel.border = element_rect(fill = NA, 
                                    size = 1)) +
  theme(strip.background = 
          element_rect(fill = "white",
                       color = "white",
                       size = 1)) +
  theme(text = element_text(size=13)) +
  theme_few()+
  xlab("Effort") + 
  ylab("CPUE") +
  ggtitle("variabless por Numero de Tripulantes y por zona")

Year x mes

ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~trim)

# No
# Zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few() +
  xlab("dfp") + ylab("Catch") +
  facet_grid(año~trim)

  facet_grid(~zona)
   <ggproto object: Class FacetGrid, Facet, gg>
       compute_layout: function
       draw_back: function
       draw_front: function
       draw_labels: function
       draw_panels: function
       finish_data: function
       init_scales: function
       map_data: function
       params: list
       setup_data: function
       setup_params: function
       shrink: TRUE
       train_scales: function
       vars: function
       super:  <ggproto object: Class FacetGrid, Facet, gg>
# CPUE slope varies between habitats - interaction
# Tipo carnada zona
ggplot(bitb3,  aes(x = dfp, y = (PESO_corr))) +
  geom_point(shape = 16, size = 3, alpha = 0.7) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
  theme_few()+
  xlab("Points") + ylab("Catch") +
  facet_grid(TIPO_CARNADA~zona)

# Perhaps

Genero categorias y Factorizo

bitb3 <- bitb3 %>%
  mutate(POTENCIA_CAT = cut(POTENCIA_MOTOR_EQ, 
                            breaks = c(0, 200, 300, 400, 1000),
                      labels = c("Bajo", "Medio-Bajo", "Medio-Alto", "Alto")),
         LAT_CAT = cut(lat_ctgr, 
                              breaks = c(-20, -35, -40, -47),
                      labels = c("Norte", "Centro", "Sur")),
         N_TRIPULANTES=factor(N_TRIPULANTES),
         año=factor(año),
         zona=factor(zona),
         trim=factor(trim)) %>% 
  drop_na(zona,
          LAT_CAT,
          prof_med_cat,
          POTENCIA_CAT) %>%
  filter(CPUE<600) %>% 
  mutate(logCPUE=log(CPUE))

glimpse(bitb3)
   Rows: 3,340
   Columns: 20
   $ N_TRIPULANTES     <fct> 4, 8, 8, 8, 6, 7, 6, 7, 8, 8, 6, 8, 8, 9, 9, 9, 9, 9…
   $ POTENCIA_MOTOR_EQ <dbl> 250, 380, 320, 450, 380, 420, 370, 240, 320, 360, 45…
   $ año               <fct> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
   $ mes               <dbl> 7, 7, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11,…
   $ trim              <fct> 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
   $ dfp               <dbl> 13.833333, 11.741667, 15.086806, 13.729167, 18.31736…
   $ lat_ctgr          <dbl> -40.66667, -39.96667, -40.50000, -33.91667, -40.6666…
   $ lon_ctgr          <dbl> -74.50000, -74.33333, -74.50000, -72.58333, -74.6666…
   $ zona              <fct> 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 3, 3, 3, 2, 2, 2, 2, 2…
   $ PESO_corr         <dbl> 1095.0568, 1256.6812, 1422.8584, 975.6000, 2026.6464…
   $ TIPO_CARNADA      <chr> "2", "3", "3", NA, "1", "2", "2", "2", "3", "3", "3"…
   $ HORAS_REPOSO      <dbl> 1400, 1800, 1200, 4800, 1200, 1200, 2400, 1200, 1200…
   $ N_PANOS           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
   $ order             <int> 8, 10, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 2…
   $ CPUE              <dbl> 79.160730, 107.027493, 94.311443, 71.060393, 110.640…
   $ CPUE2             <dbl> 0.7821834, 0.6981562, 1.1857153, 0.2032500, 1.688872…
   $ prof_med_cat      <fct> "(1e+03,1.5e+03]", "(500,1e+03]", "(1e+03,1.5e+03]",…
   $ POTENCIA_CAT      <fct> Medio-Bajo, Medio-Alto, Medio-Alto, Alto, Medio-Alto…
   $ LAT_CAT           <fct> Norte, Centro, Norte, Sur, Norte, Norte, Norte, Nort…
   $ logCPUE           <dbl> 4.371480, 4.673086, 4.546603, 4.263530, 4.706288, 4.…

The data exploration showed:

  1. One significant outlier in catch
  2. Deviation from normality in response variable
  3. Possible departure from homogeneity
  4. Few zeros in the response variable
  5. No collinearity
  6. No imbalance
  7. Possible season x habitat interaction

Chequeo distribuciones

bcpue <- ggplot(bitb3,
                aes(CPUE, group=año))+
  coord_flip()+
  geom_boxplot()+
  theme_few()


hcpue <- ggplot(bitb3 ,
                aes(CPUE))+
  geom_histogram(bins=15, fill=2)+
  theme_few()

ggarrange(bcpue, hcpue, ncol = 2)

Grafico para comprobar distribución Gamma

dist <- descdist(bitb3$CPUE, boot = 200)

Gamma distribution es la mas adecuada para la variable y que tiene la forma de :

f(x, u, \sigma) = \frac{1}{\sqrt{2hr}}e^{\frac{1(x-u)^2}{\sigma^2}}

Donde;

x \in R

En la transformación de la dsata se comprueba la inconcistencia del uso de datos de la variable CPUE transformada

Modelos GLM

#Predictores (Factores Principales)
M01 = glm(formula = CPUE ~  año,
              family = gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M02 = glm(formula= CPUE ~  año+trim,
              family = gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M03 = glm(formula= CPUE ~  año+trim+prof_med_cat,
              family =  gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M04 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES,
              family = gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M05 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT,
              family =  gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M06 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+TIPO_CARNADA,
              family =  gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M07 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona,
              family =  gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
M08 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona:prof_med_cat,
              family = gaussian(link = "identity"), 
              data = bitb3,na.action=na.exclude)
rmodelo01 <- c(AIC(M01),(M01$null.deviance-M01$deviance)/M01$null.deviance)
rmodelo02 <- c(AIC(M02),(M02$null.deviance-M02$deviance)/M02$null.deviance)
rmodelo03 <- c(AIC(M03),(M03$null.deviance-M03$deviance)/M03$null.deviance)
rmodelo04 <- c(AIC(M04),(M04$null.deviance-M04$deviance)/M04$null.deviance)
rmodelo05 <- c(AIC(M05),(M05$null.deviance-M05$deviance)/M05$null.deviance)
rmodelo06 <- c(AIC(M06),(M06$null.deviance-M06$deviance)/M06$null.deviance)
rmodelo07 <- c(AIC(M07),(M07$null.deviance-M07$deviance)/M07$null.deviance)
rmodelo08 <- c(AIC(M08),(M08$null.deviance-M08$deviance)/M08$null.deviance)
resultados <- as.data.frame(rbind(rmodelo01,rmodelo02,rmodelo03,rmodelo04,rmodelo05,rmodelo06,
                   rmodelo07,rmodelo08))

resultados <- resultados %>% 
  rename("AIC"=V1,
         "Deviance"=V2)

kbl(resultados, booktabs = T,format = "html",
    caption = "Parámetros de los modelos estimados") %>%
    kable_styling(latex_options = c("striped",
                                  "condensed","scale_down"),
                full_width = FALSE) 
predict.3    <- predict(M07,se.fit=T)
predict.mean3 <-exp(predict.3$fit)
predict.var3  <-predict.mean3^2*predict.3$se.fit^2
predict.sd3   <-sqrt(predict.var3)
predict.cv3   <-predict.sd3/predict.mean3
predict.linf3 <-predict.mean3-1.96*predict.sd3
predict.lsup3 <-predict.mean3+1.96*predict.sd3
hist(exp(predM06$fit))

Seleccion del modelo final M06

check_model(M07)

Table with significance level (***).

tab_model(M01, 
          M02,
          M03, 
          M04,
          M05,
          M06,
          M07,
          M08,
          p.style = "stars")
  CPUE CPUE CPUE CPUE CPUE CPUE CPUE CPUE
Predictors Estimates CI Estimates CI Estimates CI Estimates CI Estimates CI Estimates CI Estimates CI Estimates CI
(Intercept) 45.27 *** 22.46 – 68.08 42.90 *** 19.14 – 66.66 33.63 ** 10.36 – 56.90 -3.00 -60.39 – 54.39 31.30 -25.96 – 88.56 121.42 -44.67 – 287.50 30.55 -27.52 – 88.61 -4.34 -67.24 – 58.55
año [2014] -19.05 -44.15 – 6.04 -24.71 -50.17 – 0.76 -23.76 -48.36 – 0.83 -17.55 -42.00 – 6.90 -1.95 -26.33 – 22.43 -45.62 * -87.06 – -4.19 -3.96 -28.43 – 20.51 -8.07 -32.61 – 16.47
año [2015] -4.10 -29.54 – 21.34 -4.61 -30.14 – 20.92 -17.57 -42.24 – 7.10 -0.00 -24.87 – 24.87 14.90 -9.94 – 39.73 -25.45 -66.99 – 16.09 10.86 -14.18 – 35.90 6.08 -19.01 – 31.17
año [2016] 31.01 * 5.22 – 56.80 22.00 -4.86 – 48.87 4.41 -21.56 – 30.39 11.85 -13.89 – 37.60 26.23 * 0.63 – 51.83 -18.39 -60.45 – 23.67 21.21 -4.63 – 47.05 18.48 -7.27 – 44.23
año [2017] 6.25 -17.66 – 30.17 6.30 -17.68 – 30.28 -4.39 -27.54 – 18.77 -0.40 -23.49 – 22.68 10.37 -12.54 – 33.27 -37.39 -77.88 – 3.11 6.12 -17.10 – 29.34 2.09 -21.17 – 25.35
año [2018] -3.95 -27.68 – 19.78 -4.09 -27.86 – 19.69 -23.24 * -46.28 – -0.20 -13.79 -36.88 – 9.30 -2.77 -25.70 – 20.16 -53.80 ** -94.29 – -13.32 -10.32 -33.51 – 12.87 -15.09 -38.32 – 8.13
año [2019] 21.67 -2.77 – 46.12 21.45 -2.98 – 45.88 9.17 -14.42 – 32.77 17.39 -6.40 – 41.19 25.62 * 2.05 – 49.19 -25.36 -66.17 – 15.46 20.48 -3.08 – 44.05 14.89 -8.80 – 38.58
año [2020] -21.96 -46.10 – 2.19 -20.98 -45.22 – 3.25 -29.77 * -53.14 – -6.40 -11.82 -35.41 – 11.77 3.14 -20.37 – 26.65 -48.70 * -89.56 – -7.83 -2.38 -26.16 – 21.40 -5.47 -29.31 – 18.37
año [2021] 55.29 *** 31.14 – 79.43 49.12 *** 24.71 – 73.52 40.88 *** 17.34 – 64.42 48.56 *** 24.94 – 72.18 58.81 *** 34.95 – 82.68 8.87 -31.94 – 49.67 48.70 *** 24.62 – 72.78 48.36 *** 24.23 – 72.49
año [2022] 33.32 ** 8.43 – 58.21 27.27 * 1.95 – 52.58 20.30 -4.10 – 44.71 29.96 * 5.51 – 54.40 40.26 ** 15.55 – 64.97 -1.15 -42.46 – 40.16 31.55 * 6.67 – 56.44 29.76 * 4.87 – 54.66
trim [2] 15.86 ** 6.07 – 25.65 8.40 -1.08 – 17.88 11.33 * 1.79 – 20.87 11.17 * 1.75 – 20.60 11.09 * 1.76 – 20.43 11.02 * 1.63 – 20.41 10.47 * 1.12 – 19.82
trim [3] 19.61 ** 7.44 – 31.78 16.79 ** 5.06 – 28.53 21.93 *** 10.19 – 33.67 18.50 ** 6.89 – 30.12 14.22 * 2.70 – 25.74 17.35 ** 5.77 – 28.93 17.10 ** 5.51 – 28.70
trim [4] 1.31 -5.71 – 8.34 -2.33 -9.12 – 4.46 -1.17 -7.96 – 5.62 -1.56 -8.26 – 5.14 -4.48 -11.12 – 2.17 -0.57 -7.27 – 6.12 0.45 -6.26 – 7.17
prof med cat
[1001-1.5e+03]		 17.80 *** 11.16 – 24.44 22.47 *** 15.72 – 29.22 21.81 *** 15.14 – 28.48 23.82 *** 17.15 – 30.50 20.27 *** 13.60 – 26.95 82.41 *** 55.05 – 109.78
prof med cat [1501-2e+03] 63.66 *** 55.44 – 71.88 63.17 *** 54.86 – 71.49 59.62 *** 51.33 – 67.91 60.37 *** 52.14 – 68.61 55.48 *** 47.11 – 63.85 74.71 *** 47.06 – 102.36
prof med cat [2001-3e+03] 66.83 *** 47.76 – 85.90 68.40 *** 49.52 – 87.27 65.25 *** 46.60 – 83.90 65.75 *** 47.33 – 84.17 59.33 *** 40.65 – 78.01 64.31 ** 21.80 – 106.82
N TRIPULANTES [4] 13.90 -41.73 – 69.54 14.88 -40.06 – 69.82 17.14 -37.06 – 71.34 15.78 -38.89 – 70.45 15.63 -38.74 – 70.01
N TRIPULANTES [5] 60.36 * 6.97 – 113.75 54.09 * 1.34 – 106.85 58.17 * 6.17 – 110.18 59.18 * 6.63 – 111.72 58.84 * 6.52 – 111.16
N TRIPULANTES [6] 26.94 -25.53 – 79.42 21.62 -30.31 – 73.54 29.20 -22.01 – 80.41 38.35 -13.76 – 90.47 40.92 -11.10 – 92.93
N TRIPULANTES [7] 5.36 -47.15 – 57.86 -4.03 -56.03 – 47.96 4.67 -46.62 – 55.96 17.25 -35.16 – 69.66 19.40 -32.88 – 71.68
N TRIPULANTES [8] 31.78 -20.95 – 84.51 12.54 -39.88 – 64.96 16.01 -35.68 – 67.69 31.99 -20.74 – 84.72 33.38 -19.22 – 85.98
N TRIPULANTES [9] 37.40 -15.52 – 90.33 3.70 -49.26 – 56.66 14.32 -37.94 – 66.58 21.57 -31.54 – 74.68 23.66 -29.29 – 76.60
N TRIPULANTES [10] 1.39 -56.22 – 59.00 -39.17 -96.77 – 18.43 -27.89 -84.70 – 28.91 -24.82 -82.69 – 33.05 -22.19 -80.19 – 35.81
LAT CAT [Centro] -31.00 *** -40.16 – -21.83 -26.66 *** -35.84 – -17.48 -9.65 -25.97 – 6.66 -8.56 -24.85 – 7.72
LAT CAT [Sur] -43.80 *** -53.01 – -34.59 -37.90 *** -47.13 – -28.66 -29.35 *** -46.16 – -12.55 -26.63 ** -43.40 – -9.85
TIPO CARNADA [1] -87.20 -239.67 – 65.28
TIPO CARNADA [2] -50.33 -202.64 – 101.98
TIPO CARNADA [3] 59.29 -95.96 – 214.55
zona [2] -30.26 *** -42.07 – -18.44
zona [3] -6.41 -26.55 – 13.72
prof med cat [501-1e+03]
× zona2		 6.93 -18.62 – 32.47
prof med cat
[1001-1.5e+03] × zona2		 -62.36 *** -77.70 – -47.01
prof med cat [1501-2e+03]
× zona2		 -10.01 -27.22 – 7.20
prof med cat [2001-3e+03]
× zona2		 -0.40 -43.42 – 42.62
prof med cat [501-1e+03]
× zona3		 24.17 -7.32 – 55.66
prof med cat
[1001-1.5e+03] × zona3		 -32.27 ** -55.94 – -8.60
prof med cat [1501-2e+03]
× zona3		 13.12 -12.43 – 38.67
prof med cat [2001-3e+03]
× zona3		 61.51 * 6.00 – 117.02
Observations 3340 3340 3340 3340 3340 3293 3340 3340
R2 0.078 0.083 0.150 0.176 0.197 0.229 0.205 0.216
* p<0.05   ** p<0.01   *** p<0.001

Table comparing performance model.

compare_performance(M01, 
          M02,
          M03, 
          M04,
          M05,
          M06,
          M07,
          M08,
          rank = TRUE, 
          verbose = FALSE)
   # Comparison of Model Performance Indices
   
   Name | Model |    R2 |   RMSE |  Sigma | Performance-Score
   ----------------------------------------------------------
   M06  |   glm | 0.229 | 77.035 | 77.364 |           100.00%
   M08  |   glm | 0.216 | 77.282 | 77.667 |            94.44%
   M07  |   glm | 0.205 | 77.804 | 78.120 |            87.21%
   M05  |   glm | 0.197 | 78.206 | 78.500 |            81.47%
   M04  |   glm | 0.176 | 79.230 | 79.504 |            66.63%
   M03  |   glm | 0.150 | 80.476 | 80.670 |            48.76%
   M02  |   glm | 0.083 | 83.567 | 83.731 |             3.21%
   M01  |   glm | 0.078 | 83.796 | 83.921 |             0.00%

Plot comparing performance model.

plot(compare_performance(M01, 
          M02,
          M03, 
          M04,
          M05,
          M06,
          M07,
          M08,
          rank = TRUE, 
          verbose = FALSE))
model_performance(M06)
   # Indices of model performance
   
   R2    |   RMSE |  Sigma
   -----------------------
   0.229 | 77.035 | 77.364
res_glm<-as.numeric(M06$residuals)

res_hist<-ggplot(as.data.frame(res_glm), aes(x=res_glm)) +
  geom_histogram(position="identity", bins = 20)+
  theme_bw()+
  xlab("Residuos")+
  ylab("Frecuencia")
#res_qqnorm<-qqnorm(res_glm)+
 # qqline(res_glm) 

res_density<-ggdensity(res_glm, res_glm = "", fill = "lightgray", 
                                title = "") +
  scale_x_continuous() +
  stat_overlay_normal_density(color = "red", linetype = "dashed")+
  xlab("Residuos")+
  ylab("Densidad")
res_points<-ggplot(as.data.frame(res_glm), 
                            aes(y=res_glm, x=(1:length(res_glm)))) +
  geom_point()+
  xlab("")+
  ylab("Residuos")+
  theme_few()
ggarrange(res_hist,  res_density, res_points,
          labels = c("A", "B", "C"),
          ncol = 3)

output<-M01$fitted.values

Regarding this result, we can select best model and check some assumtion like:

  • Colinealidad
  • Data points influyentes
  • Homoscedasticidad
  • Normalidad de los residuos
  • Independencia

COMPOSICIONES DE TALLAS

La información contenida en las estructiras de tallas son consideradas una de las piezas mas importantes en la ciencia pesquera (Canales, Punt, & Mardones, 2021; Hordyk, Ono, Sainsbury, Loneragan, & Prince, 2014; Rudd & Thorson, 2018). Sin embargo, esta debe ser validada y confirmada para no violar supuestos referidos a la repreentatividad del dato y su relación con la dinámica poblacional.

Ahora nos disponemos a explorar esta fuente de información.

Primero identificamos la estructura de la base;

glimpse(long)
   Rows: 15,831
   Columns: 24
   $ ID                  <chr> "COD_BARCO=730881 AND TO_CHAR(FECHA_HORA_RECALADA,…
   $ COD_BARCO           <dbl> 730881, 730881, 730881, 730881, 730881, 730881, 73…
   $ FECHA_HORA_RECALADA <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
   $ FECHA_HORA_ZARPE    <dttm> 2021-05-10 11:00:00, 2021-05-10 11:00:00, 2021-05…
   $ COD_PESQUERIA       <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50…
   $ PUERTO_RECALADA     <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 3,…
   $ NRO_FORMULARIO      <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,…
   $ REGION              <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1,…
   $ FECHA_LANCE         <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
   $ NUMERO_LANCE_EX     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
   $ LATITUD             <dbl> 243900, 243900, 243900, 243900, 243900, 243900, 24…
   $ LONGITUD            <dbl> 714400, 714400, 714400, 714400, 714400, 714400, 71…
   $ NUMERO_EJEMPLARES   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ PESO_TOTAL_MUESTRA  <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ NUMERO_CAJA         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
   $ COD_ESPECIE         <dbl> 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37…
   $ ORIGEN_MUESTRA      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3,…
   $ FECHA_MUESTREO      <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ N_TOTAL_INDIV       <dbl> 85, 85, 85, 85, 85, 85, 85, 84, 84, 84, 84, 84, 23…
   $ LONGITUD_DESCARTE   <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
   $ LONGITUD_MUESTRA    <dbl> 115, 125, 129, 130, 143, 152, 157, 100, 108, 110, …
   $ SEXO                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
   $ N_INDIVIDUOS        <dbl> 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 1, 1,…
   $ año                 <dbl> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 20…
colSums(is.na(long))
                    ID           COD_BARCO FECHA_HORA_RECALADA    FECHA_HORA_ZARPE 
                  2997                   0                   0                   0 
         COD_PESQUERIA     PUERTO_RECALADA      NRO_FORMULARIO              REGION 
                     0                   0                   0                   0 
           FECHA_LANCE     NUMERO_LANCE_EX             LATITUD            LONGITUD 
                     0                   0                4425                4425 
     NUMERO_EJEMPLARES  PESO_TOTAL_MUESTRA         NUMERO_CAJA         COD_ESPECIE 
                 15733               15831                   0                   0 
        ORIGEN_MUESTRA      FECHA_MUESTREO       N_TOTAL_INDIV   LONGITUD_DESCARTE 
                     0               15831                   0               15831 
      LONGITUD_MUESTRA                SEXO        N_INDIVIDUOS                 año 
                     0                   0                   2                   0

Primero selecciono las columnas de interes y luego genero la expansión de LONGITUD_MUESTRA a N_INDIVIDUOS. Expand frecuency data related length, in this case N_INDIVIDUOS column have frecuency that we need expand to whole data frame.

Selecciono variables de interes

long1 <- long %>% 
  dplyr::select(6,8,11,12,21,22,23,24)
dim(long1)
   [1] 15831     8
long2 <- long %>% 
  drop_na(N_INDIVIDUOS) %>% 
  type.convert(as.is = TRUE) %>% 
  uncount(N_INDIVIDUOS)

dim(long2)
   [1] 34299    23
table(long2$REGION)
   
       1     3     4     5     7     8    10    12    14 
     494  3055  1697 12111  9452  4263  2848     1   378

legend.labels <- c('Indet.', 'Macho' , 'Hembra')
nbco <- ggplot(long2 %>% 
                 filter(REGION!=12), aes(x=LONGITUD_MUESTRA, y = as.factor(año),
                         fill=as.factor(SEXO)))+
  geom_density_ridges(stat = "density_ridges", bins = 20, 
                      scale = 2, draw_baseline = FALSE,
                      alpha=0.5)+
  facet_wrap(.~REGION, ncol=4) +
  geom_vline(xintercept = 110, color = "red")+
  scale_fill_viridis_d(option = "G",
                       name="SEXO",
                       labels=legend.labels)+
  scale_y_discrete(breaks = seq(from = 2004, to = 2022, by = 2))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  theme_few()+
  xlab("Longitud (cm.)")+
  ylab("")
#scale_x_discrete((limits = rev(levels(talla2021$ANO_ARR))))+
nbco

Prepara los vectores para sumar a los .dat del modelo si los necesito.

long2$LONG_CAT <- as.numeric(as.character(cut(x = long2$LONGITUD_MUESTRA, 
                                              breaks = seq(10,220,2), 
                                              labels = seq(10,218,2),
                                              right = FALSE)))
LONGT <- table(long2$año, long2$LONG_CAT)
# guardar en un .xls

BIOLOGICO

names(bio)
    [1] "ID"                             "COD_BARCO"                     
    [3] "FECHA_HORA_RECALADA"            "FECHA_HORA_ZARPE"              
    [5] "PUERTO_RECALADA"                "PUERTO_ZARPE"                  
    [7] "REGION"                         "COD_PESQUERIA"                 
    [9] "NRO_FORMULARIO"                 "FECHA_LANCE"                   
   [11] "NUMERO_LANCE_EX"                "LATITUD"                       
   [13] "LONGITUD"                       "NUMERO_EJEMPLARES"             
   [15] "PESO_TOTAL_MUESTRA"             "COD_ESPECIE"                   
   [17] "N_ESPECIMEN"                    "ORIGEN_MUESTRA"                
   [19] "SEXO_ESPECIMEN"                 "NRO_MUESTRA"                   
   [21] "LONGITUD_ESPECIMEN"             "PESO_ESPECIMEN"                
   [23] "NUMERO_OTOLITOS"                "PESO_EVISCERADO"               
   [25] "MADUREZ"                        "PESO_GONADAS"                  
   [27] "CONTENIDO_ESTOMACAL"            "REPLECION"                     
   [29] "LONGITUD_DISCO"                 "ANCHO_DISCO"                   
   [31] "PESO_ESTOMAGO"                  "PESO_PARED_ESTOMAGO"           
   [33] "NUMERO_FICHA"                   "CONT_ESTOMACAL_ESTADO_ESTOMAGO"
   [35] "NUMERO_FICHA_GONADA"            "año"
table(bio$COD_ESPECIE)
   
       37 
   127409

Filtro ciertos datos para estimar coeficientes

bio2 <- bio %>% 
  drop_na(PESO_ESPECIMEN,
          LONGITUD_ESPECIMEN,
          SEXO_ESPECIMEN) %>% 
  filter(PESO_ESPECIMEN>100,
         año>2005)

Ahora por sexo

legend.label <- c('Macho' , 'Hembra')
my.plotsex<- ggplot(bio2 %>% 
                      filter(SEXO_ESPECIMEN!="0",
                             SEXO_ESPECIMEN!="3",
                             SEXO_ESPECIMEN!="6"), 
                  aes(y=PESO_ESPECIMEN, x=LONGITUD_ESPECIMEN, 
                      colour=SEXO_ESPECIMEN)) +
  geom_point(alpha=.1) + 
  scale_colour_manual(values = c("black", "red"),
                        name="Sexo",
                      labels=legend.label)+
  facet_wrap(~ año, ncol = 4) + 
  geom_smooth(method = "loess", 
              se=TRUE,  
              formula = y ~ x, 
              size=0.3,
              span=5)+
  theme_few()+
  theme(legend.position = "bottom")+
  labs(y="Peso (gr.)",
       x="Longitud (cm)")
my.plotsex

Ahora por región

legend.label <- c('Macho' , 'Hembra')
my.plotreg<- ggplot(bio2 %>% 
                      drop_na(REGION) %>% 
                      filter(SEXO_ESPECIMEN!="0",
                             SEXO_ESPECIMEN!="3",
                             SEXO_ESPECIMEN!="6",
                             REGION!=12), 
                  aes(y=PESO_ESPECIMEN, x=LONGITUD_ESPECIMEN, 
                      colour=SEXO_ESPECIMEN)) +
  geom_point(alpha=.1) + 
  scale_colour_manual(values = c("black", "red"),
                        name="Sexo",
                      labels=legend.label)+
  facet_wrap(~ REGION, ncol = 4) + 
  geom_smooth(method = "loess", 
              se=TRUE,  
              formula = y ~ x, 
              size=0.3,
              span=5)+
  theme_few()+
  theme(legend.position = "bottom")+
  labs(y="Peso (gr.)",
       x="Longitud (cm)")
my.plotreg

Extraigo coeficientes globales

x<-log(bio2$LONGITUD_ESPECIMEN)
y<-log(bio2$PESO_ESPECIMEN)
reg.l<-lm(y~x)
summary(reg.l)
   
   Call:
   lm(formula = y ~ x)
   
   Residuals:
        Min       1Q   Median       3Q      Max 
   -2.68627 -0.07255 -0.00070  0.07109  2.47970 
   
   Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
   (Intercept) -5.005673   0.018393  -272.1   <2e-16 ***
   x            3.082291   0.004136   745.3   <2e-16 ***
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   
   Residual standard error: 0.1557 on 31985 degrees of freedom
   Multiple R-squared:  0.9456, Adjusted R-squared:  0.9455 
   F-statistic: 5.555e+05 on 1 and 31985 DF,  p-value: < 2.2e-16
anova(reg.l)
   Analysis of Variance Table
   
   Response: y
                Df  Sum Sq Mean Sq F value    Pr(>F)    
   x             1 13468.9   13469  555453 < 2.2e-16 ***
   Residuals 31985   775.6       0                      
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res <- residuals(reg.l)
plot(res)
abline(h = 0, col = "red", lty=2)

BITACORA FIPA 96- 14


bitfipa<- bit8696 %>% 
  mutate(CPUE = CAPTURA/DFP) %>% 
  drop_na(CPUE)


cpuefipa <- ggplot(bitfipa %>%
                     filter(ANO<1997)  %>%
                  group_by(ANO) %>% 
                  summarise(CPUEF=mean(CPUE)),
                aes(ANO, CPUEF))+
  geom_point(shape = 16, size = 3, alpha = 0.7)+
  scale_x_continuous(breaks = seq(from = 1986, to = 2022, by = 1))+
  geom_smooth(method = 'lm', 
              colour = 'blue', 
              size = 1.5)+
  # stat_regline_equation(label.x=1993, label.y=500)+
  # stat_cor(aes(label=..rr.label..), label.x=1993, label.y=450)+
  #facet_wrap(~AREA)+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2),
        panel.grid = element_blank())+
  labs(y="CPUE (kg/dfp)",
       x="")
cpuefipa

Guardo datos de la CPUE nominal del FIPA 96-14

cpuenomfipa <- bitfipa %>%
  filter(ANO<1997)  %>%
  group_by(ANO) %>% 
  summarise(CPUEF=mean(CPUE))
write_csv(cpuenomfipa, "CPUE_FIPA.csv")

MAPAS

Lo primero es transformar los datos de lon_ctgr y lat_ctgr a formarto geom. A su vez, saco las NA de las coord

transformar los datos en un sf object

codmap <- st_as_sf(bitb2 %>% 
                    drop_na(lon_ctgr) %>% 
                    drop_na(lat_ctgr), 
                   coords = c("lon_ctgr", "lat_ctgr"),  crs = 4326)

Ahora genero los mapas de Chile con un raster.

chile <- raster::getData("GADM", country = "CHL", level = 0)
chile1<-fortify(chile)
chilemap <- ggplot()+
   geom_polygon(data=chile, aes(x=long, y=lat, group=group),
              fill="lightblue",color="grey20", size=0.15)+
   coord_sf(crs = st_crs(4326),
            xlim = c(-80, -65),
            ylim = c(-58, -15)) +
          theme_bw()
 chilemap

# Aca veo los nombres
chile@data$NAME_1
   NULL

Luego genero los bordes sobre los cuales haré la grilla.

e <- extent(-80,-65,-58,-15)
rc <- crop(chile, e)

# la proyección adecuada en lat long
tcrs <- CRS("+init=epsg:4326")
transformed_points <- spTransform(rc, tcrs)
# para dejarlo en formato geom_sf
chile2 <- st_as_sf(transformed_points)

Se estructura la grilla en el objeto raster de chile1

grid <- chile2 %>%
  st_make_grid(cellsize = c(1,1)) %>% # para que quede cuadrada
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>% # objeto en spatial feature
  mutate(cellid = row_number())

Pongo los datos codmap en la grilla. Aca solo elegí dfpy CPUE pero se ùeden resumir otros

joindat <- grid %>%
  st_join(codmap) %>% 
  group_by(cellid, año) %>% 
  summarise(CPUEM = mean(CPUE),
            DFPM =mean(dfp)) # por ejemplo, plotear la captura "CAPTURA_1"

Mapa Esfuerzo

here we can viz diferent variables

## Plot final
mas1 <- ggplot() +
  geom_sf(data=joindat %>% 
             filter(!is.na(DFPM)), aes(fill = DFPM),
          color=NA) +
  scale_fill_viridis_b(option="E",
                       direction=-1, name="DFP")+
  geom_sf(data = grid,  fill=NA, color=NA) +
  geom_sf(data = chile2, color="grey", fill="white") +
  coord_sf() +
  geom_hline(yintercept = -47, color = "red")+
  scale_alpha(guide="none")+
  facet_wrap(~año, ncol=6)+
  scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides(colour = guide_legend()) +
  theme_bw()+
  theme(panel.background = element_rect(fill = 'white'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())
mas1

Mapa CPUE

## Plot final
mas2 <- ggplot() +
  geom_sf(data=joindat %>% 
             filter(!is.na(CPUEM)), aes(fill = CPUEM),
          color=NA) +
  scale_fill_viridis_b(option="G",
                       direction=-1, name="CPUE")+
  geom_sf(data = grid,  fill=NA, color=NA) +
  geom_sf(data = chile2, color="grey", fill="white") +
  coord_sf() +
  geom_hline(yintercept = -47, color = "red")+
  scale_alpha(guide="none")+
  facet_wrap(~año, ncol=6)+
  scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides(colour = guide_legend()) +
  theme_bw()+
  theme(panel.background = element_rect(fill = 'white'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())
mas2

Falta componer mapas de tallas medias y alguna otra variable de interés.

En función de los datos analizados, el modelo propuesto es un enfoque de producción global estado espacio descrito por Pedersen & Berg (2017) y Acom (2015).

REFERENCIAS

Acom, I. C. M. (2015). ICES WKLIFE V REPORT 2015 Methodologies based on Life-history Traits , Development of Quantitative Assessment Exploitation Characteristics and other Report of the Fifth Workshop on the Relevant Parameters for Data-limited Stocks ( WKLIFE V ) Lisbon , Port. (October), 5–9.
Canales, C. M., Punt, A. E., & Mardones, M. (2021). Can a length-based pseudo-cohort analysis (LBPA) using multiple catch length-frequencies provide insight into population status in data-poor situations? Fisheries Research, 234(October 2020), 105810. https://doi.org/10.1016/j.fishres.2020.105810
Hordyk, A., Ono, K., Sainsbury, K., Loneragan, N., & Prince, J. (2014). Some explorations of the life history ratios to describe length composition, spawning-per-recruit, and the spawning potential ratio. ICES Journal of Marine Science, 72(1), 204–216. https://doi.org/10.1093/icesjms/fst235
Pedersen, M. W., & Berg, C. W. (2017). A stochastic surplus production model in continuous time. Fish and Fisheries, 18(2), 226–243. https://doi.org/10.1111/faf.12174
Rudd, M. B., & Thorson, J. T. (2018). Accounting for variable recruitment and fishing mortality in length-based stock assessments for data-limited fisheries. Canadian Journal of Fisheries and Aquatic Sciences, 75(7), 1019–1035. https://doi.org/10.1139/cjfas-2017-0143